arXiv: 1508.02600V1 [math.NA] 30Jul2015 


An adaptive multiresolution method for ideal 
magnetohydrodynamics using divergence cleaning with 
parabolic-hyperbolic correction 


Anna Karina Fontes Gomes'^’^’®, Margarete Oliveira Domingnes^’'^’*^, 
Kai Schneiderf Odim Mendes'^’^’®, Ralf Deiterding® 

°‘P6s-Graduagao em Computagdo Aplicada (CAP) 

^Laboratorio Associado de Computagdo e Matemdtica Aplicada (LAC) 
‘^Coordenadoria dos Laboratorios Associados (CTE) 

‘^Divisdo de Ceofi'sica Espacial, Coordenagdo de Ciencias Espaciais(CEA) 
^Instituto Nacional de Pesquisas Espaciais (INPE), Av. dos Astronautas 1758, 
12227-010 Sdo Jose dos Campos, Sdo Paulo, Brazil 
^M2P2-CNRS & Centre de Mathematiques et d’Informatique (CMI), Aix-Marseille 
Universite, 38 rue F. Joliot-Curie, 13)51 Marseille Cedex 20, France 
^German Aerospace Center (DLR), Institute of Aerodynamics and Flow Technology, 
Bunsenstr. 10, 37073 Gottingen, Germany 


Abstract 

We present an adaptive multiresolntion method for the nnmerical simulation 
of ideal magnetohydrodynamics in two space dimensions. The discretiza¬ 
tion uses a hnite volume scheme based on a Cartesian mesh and an explicit 
compact Runge-Kutta scheme for time integration. Harten’s cell average 
multiresolution allows to introduce a locally refined spatial mesh while con¬ 
trolling the error. The incompressibility of the magnetic field is controlled 
by using a Generalized Lagrangian Multiplier (GLM) approach with a mixed 
hyperbolic-parabolic correction. Different applications to two-dimensional 
problems illustrate the properties of the method. For each application CPU 
time and memory savings are reported and numerical aspects of the method 
are discussed. The accuracy of the adaptive computations is assessed by 
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comparison with reference solutions computed on a regular fine mesh. 

Keywords: Magnetohydrodynamics, Multiresolution Analysis, Finite 
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1. Introduction 


The magnetohydrodynamic (MHD) equations, which consist of the com¬ 
pressible Euler equations of hydrodynamics coupled with the Maxwell equa¬ 
tions of electrodynamics, are used for mathematical modeling of numerous 
phenomena encountered in our daily life. Prominent examples can be found 
in the physics of the Sun-Earth’s electrodynamical interaction chain, and in 
the dynamo action caused by motion of liquid metal inside the mantle of 
the Earth, which generates its magnetic held. The numerical challenge for 
solving the ideal MHD equations, a coupled set of nonlinear Partial Differ¬ 
ential Equations (PDEs), is the presence of multiple spatial and temporal 
scales. The complex character of boundary conditions of the magnetic held, 
in comparison to that one for the classical hydrodynamics, requires even more 
sophisticated approaches. In a surrounding vacuum, for example, the mag¬ 
netic held does not vanish, it only decays. Thus, at the boundary it has to be 
matched with the held of the huid region. A second difficulty is to maintain 
the incompressibility of the magnetic held numerically, which is imposed by 
Gauss’ law. Therefore, in the numerical simulations, special attention has to 
be paid to this incompressibility, because, as shown in practice, uncontrolled 
divergence errors can modify the underlying physics. For details we refer 
the reader to, e.g., 0,H 1^. [^ . Typically, projection methods based on 
the Helmholtz decomposition are used. These methods are computationally 
demanding, especially in three-dimension, because the solution of an elliptic 
problem requires a Poisson equation solver. An alternative method is the 
divergence cleaning one, which is based on Lagrangian multipliers. In the 
hnite element context, Assous et ah [ij introduced this approach for time- 
dependent Maxwell equations. Several variants can be found in the literature 
0 , 1 , 00 - 

In the current paper we apply the multiresolution approach to an ideal 
MHD numerical model called the Generalized Lagrange Multiplier (GLI^ 
with a mixed hyperbolic-parabolic correction proposed by Dedner et ah j6| 
to deal with the magnetic held incompressibility condition. The ideas of the 
Lagrangian multiplier formulation in this context were introduced by Munz et 
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al. [2^ in the context of Maxwell equations. With the motivation to reduce 


CPU time and memory requirements, we use an auto-adaptive discretiza¬ 
tion which is based on the multiresolution representation. The underlying 
time dependent conservation laws are discretized with hnite volume schemes 
and local grid rehnement is triggered by multiresolution analysis of the cell 
averages and thresholding of the resulting coefficients. The adaptive rehne¬ 
ment /mesh tracks steep gradients in the solution of the equation and allows 
automatic error control. For reviews on multiresolution techniques for PDFs 
we refer to [13, 18, [^, 131 and references therein. 

Preliminary results for a quasi-one dimensional MHD Riemann problem 
with exact solution have been presented in 1^, which showed the feasibil¬ 
ity of using adaptive discretizations and magnetic held divergence cleaning 
for extended GLM-MHD with local and controlled time methods. In its 
extended form, source terms similar to those in are introduced. The 
starting point is the adaptive multiresolution code originally developed by 
Roussel et al. |27| in which the Maxwell equations governing the magnetic 


held have been included [l^. In the present work, we have chosen the GLM- 
MHD approach instead of its extended version, because the divergence errors 
and the solution obtained for both cases are almost the same for the studied 
problem. A similar choice is suggested in the conclusion in [^ . The resulting 
new method has been applied to a two-dimensional Riemann test problem, 
for which a reference solution on a hne grid has been computed. The ac¬ 
curacy of the adaptive computations has been assessed and their efficiency 
in terms of memory compression compared to a hnite volume scheme on a 
regular grid has been analyzed. 

The paper is organized as follows: After a presentation of the governing 
ideal MHD equations in Section HI we recall the divergence cleaning tech¬ 
nique based on the GLM formulation in Section |2i In Section 0] space and 
time discretizations are briehy described together with the GLM discretiza¬ 
tion. In Section |5l numerical results are presented. In the last section, some 
conclusions are drawn and perspectives for future work are presented. 


2. Governing equations 

The ideal magnetohydrodynamics equations describe the dynamics of a 
compressible, inviscid and perfectly electrically conducting huid interacting 
with a magnetic held, see, e.g. [l^. The equations combine the Euler equa¬ 
tions with the Maxwell equations. The latter yields an evolution equation for 
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the magnetic field, also called induction equation, and an incompressibility 
constraint using Gauss’ law. The system of MHD equations is given by 

dp 

— + V • (pu) = 0, (Mass conservation) (la) 


dE _ 


B • B 

E + p-\ -— ) u — (u • B) B 


dpu. 


+ v 


pu*u + (p + 


B B 


I-B*B 


= 0, (Energy conservation) (lb) 
= 0, (Momentum conservation) (Ic) 


dB 


+ V • (u*B — B*u) = 0, (Induction equation) (Id) 


where p represents density, p the pressure, u = {ux,Uy,Uz) the velocity vec¬ 
tor, B = {Brc, By, Bz) the magnetic field vector, and t denotes the transposi¬ 
tion. The identity tensor of order 2 is denoted by I (the unit dyadid, that here 
corresponds to the unit matrix 3x3), and 7 the adiabatic constant (77 1 ). 
The pressure is given by the constitutive law p = (7 — 1) (£' — . 

The above system is completed by suitable initial and boundary conditions. 
In this paper this system is considered in its two-dimensional form, i.e., the 
quantities depend on two variables only {x and y). 

In this classical MHD model, the magnetic field has to satisfy the diver¬ 
gence constraint 

V ■ B = 0. (2) 

which implies the non-existence of magnetic monopoles. By rewriting the 

5B 

induction equation, we have -|- V x (B x u) = 0. Therefore, the appli¬ 
ed 

cation of the divergence operator yields — (V ■ B) = 0, as V ■ (Vx ) = 0. 
This formulation shows that if the initial condition of the magnetic field is 
divergence-free, the system will remain divergence-free along the evolution. 
However, numerically the incompressibility of the magnetic field is not nec¬ 
essarily preserved, and thus, non-physical results could be obtained or the 
computations may even become unstable Q. Since the 1980ies typical nu¬ 
merical MHD methodologies consider the enforcement of the divergence-free 
constraint. There are many techniques to perform the divergence cleaning 
in the MHD numerical models j^. In the context of this study, we have in 
mind the application of the multiresolution method based on a finite volume 
discretization with explicit time integration. Thus, the technique developed 
in Dedner et al.j^ called GLM-MHD with the mixed parabolic-hyperbolic 
correction, is well suited. Details are given in the next section. 
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3. Generalized Lagrangian multipliers for divergence cleaning 


Dedner et al. [6| proposed the GLM formulation with the hyperbolic- 
parabolic correction. Its implementation into a pre-existing MHD model is 
straightforward. An additional scalar held ijj is introduced, which couples the 
divergence constraint equation (Eq. [2]) to Faraday’s law, modifying the in¬ 
duction equation fEq. IldD . Moreover, some source terms are added similarly 
to what was proposed in The model contains one parameter related 

to the hyperbolic correction, namely c^, responsible for the propagation of 
the divergence errors, and another one related to the parabolic correction Cp, 
responsible for the damping of the monopoles. The remaining terms in the 
equations remain unchanged. The conservative characteristic of this system 
is not lost for the GLM approach. 

The resulting GLM-MHD equations written in two-dimensional form read 


dp ^ dpux , dpuy 
dy 


+ 


dt dx 
dE d 
dt dx 
dp 
dy 
d (pug;) 

dt dx 


= 0 , 


B • B' 

E + p^ -^- 1 Ux - (u • B) Brr 


E + P + 


2 

B B 




Uy-(u-B)Bj^ 


B B 


= 0 , 

d 


(3a) 


(3b) 


pul+p[p^ -^ ]-BI {pUxUy - B^By) = 0,(3c) 


d{pUy) d r, \ d 2 

+ ^ipUa,Uy-B^By) + — pUy+p[p + 


dt 


B B 


-Bf. 


d (puz') d d 

~dx ^zBx) ~\~ {,P'^z'^y BzBy) — 0, 

dBx d'lb d . ^ ^ ^ 

dBy d. n ^ r> 

^ ^ i^-By - B^Uy) + ^ - 0> 

dB d 3 

^ ~3x BzUx) “h (^VjyBz Byllz) — 0, 

dt ^ \ dx dy ) Cp ’ 


= 0,(3d) 
(3e) 
(3f) 
(3g) 
(3h) 
(3i) 


where B ■ B = u ■ B — ti^Bx ~\~ 'UyBy Cp cLnd. Ch RrG 
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the parabolic-hyperbolic parameters, with > 0. In j6|] it is defined as 


Ch = Ch(t) := CcFL 


min{Ax, Ay} 
At 


where ccfl G (0,1), Ax and Ay are the space step in x— and ?/—direction, 
respectively. At is the time step. If the parameter Ch is defined, as for instance 
in Eq. |3l then Cp is a free parameter in Eq. [3il We follow a choice proposed 
in to avoid that Cp is strongly dependent on the mesh size and the scheme 
used. Their numerical experiments showed that choosing c^/ch = 0.18, mir¬ 
rors properly the ratio between hyperbolic and parabolic effects. With this 
choice in the one-dimensional case the damping of the divergence errors oc¬ 
curs on the time scale Cp\/t and the transport of the divergence errors to the 
boundary takes place on the time scale Cht (as discussed in j^. Appendix 
A. 16 and A.l^. However, other possible choices of these parameters can be 
found in j^, and for the CTU-GLM approach in . 

Considering the vector of conservative quantities Q = (p, E, pu, B, ■0), 
the GLM-MHD system could be written compactly as 

^ + V-F(Q) = S(Q), 

where F(Q) is the physical flux, and S(Q) contains all source terms. 


4. Adaptive space and time discretization 

A finite volume discretization of the GLM-MHD system is applied, which 
results in a system of ordinary differential equations. Approximate solutions 
at a sequence of time instants are obtained by using an explicit ordinary 
differential equation solver. Here, an explicit Runge-Kutta scheme of second 
order is used. 

In the GLM-MHD Finite Volume (FV) reference scheme, we consider 
the initial value of the variable xfj as zero. The parameter Ch has a strong 
influence in the correction. In each time step, we compute the parameter 
Ch, then the GLM-MHD system is solved. First, a dimensional splitting is 
performed in x-direction, where the fluxes in the interface are treated and 
the solution updated. This procedure follows the steps: 

1. The component of the magnetic field in the x-direction flux fEq. l3fj) . 
and the divergence constraint equation (Eq. |3l]), are decoupled from the 
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other variables. These two equations form the system 


dBx dip 
dt ^ dx 


= 0 , 


dip 

'm 


4 


dB, 

dx 


rV', 


( 4 ) 

( 5 ) 


such that the local Riemann problem can be solved analytically, where 
the numerical flux in the interface is {ipm, 4i^x,m) for B^ and ip. Simi¬ 
larly as what is described in j^, we have 


f B^^ra \ ^ ( Bx,L \ , f ~ Bx,l) “ ^ii’R ~ \ 

\ J Ww V U^R-i^L)-^-tiBx,R-Bx,L) ) 

where the sub-index L, R are related to the left or right-hand state. 


2. Therefore, the numerical flux is evaluated in two steps. First we com¬ 
pute the numerical flux not considering the B^ and ip equations as 
described above, then we add the numerical flux in the interface. In 
this work, we use the Harten-Lax-van Leer-Discontinuities numerical 
flux (HLLD) with four intermediary states Qpp, and Q^, di¬ 
vided by the waves with speed Sl, Sm, e Sr, as discussed in 
the Appendix A. The states Q* and Q** are dehned as 


Qa = ipl,Ep,p] 




,,C) and Q:* = (p: 


TT’A* -.I 

Ea , Pa ^ 




with a denoting left (L) or right (i?) states. 


3. The same procedure is done for By in the ^/-direction. 


4. The computed values of ip are used to update the mixed correction 
source term for ip"‘~^^, computing -0”+^ = exp ip. 

The adaptive Multiresolution (MR) method of the present paper has been 
designed to speed up finite volume schemes for conservation laws. In the fol¬ 
lowing, a brief summary of this technique is given. For a detailed description 


of these strategies, we refer to [27|, llJ, ll2|, lUl, ll3 
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The key ingredient of MR schemes is the decay properties of the wavelet 
coefficients of the numerical solution. The decay rate indicates the local 
regularity of the solution. In regions where the solution is smooth the coeffi¬ 
cients are of small magnitude and thus coarser meshes can be used. In regions 
where the coefficients are signihcant the numerical solution is less smooth and 
strong gradients or even jumps are present and a hne mesh must be used j^. 
Stopping the rehnement in a cell at a certain scale level, where the wavelet 
coefficients are non-signihcant leads to an adaptive MR representation. 

For a hnite volume scheme the uniform cell-average representation is re¬ 
placed by cell-averages on an adaptive locally rehned mesh, which is formed 
by the cells whose wavelet coefficients are signihcant and above a given 
threshold. An example of an adaptive Cartesian mesh is presented in Fig. [T] 
In MHD solutions localized structures are present, such as discontinuities or 



Figure 1: Example of a zoom in a dyadic adaptive Cartesian mesh. Regions where the 
mesh is refined are associated with detected structures in the solution, ie., where the 
wavelet coefficients are significant. 

shocks. They could appear in different space positions in different variables. 
Thus, the adaptive mesh of the MHD system is a union of the individual 
adaptive meshes of each quantity. 

Tree structures are the natural way to store the reduced MR data. Mesh 
adaptivity is then related to an incomplete tree and the rehnement can be 
interrupted at intermediate scale levels. In other words, using the tree ter¬ 
minology, a MR mesh is formed by leaves, which are nodes without children. 
These leaves correspond to the cell which is being evolved in time. In sum¬ 
mary, there are three steps in the application of a MR scheme: rehnement, 
evolution, and coarsening. The rehnement operator accounts for possible 
translations of the solution or the creation of hner scales in the solution be¬ 
tween two subsequent time steps. Since the localized structures and thus the 









local regularity of the solution may change with time, the MR mesh at time 
t"' may not be sufficient any more at the next time step Hence, before 
evolving the solution in time, the representation of the solution should be 
interpolated onto an extended mesh that is expected to be a rehnement of 
the adaptive mesh at t"', and to contain the adaptive mesh at After 

that, the time evolution operator is applied to the leaves of the extended 
mesh. The numerical fluxes between cells of different levels are computed by 
adding extra cells, called virtual leaves, which will however not be used in 
the time evolution. Conservation is ensured by the fact that the fluxes are 
always computed on a higher level, the value being projected onto the leaves 
of a lower level. Then, wavelet thresholding is applied in order to unrehne the 
cells in the extended grid (coarsening) that are not necessary for an accurate 
representation of the solution at This data compression is based on 

the dehnition of deletable cells, where the wavelet coefficients which are not 
signihcant, he., their magnitudes are below a threshold parameter e^, where 
^ denotes the cell scale level, are called deletable cells. The data compression 
is the given by 

N 

100 

Dc = --, 

2Ln 


where N is the total number of iterations and Cn{i) is the number of cells 
in the adaptive mesh at iteration i G {I,-- - ,N}. The number of cells on 
the hnest mesh is dehned as 2^, where L the hnest scale level. However, to 
compute the flux in a conservative form, additional neighbor cells at the same 
level are also necessary. These neighbor cells are not necessarily present on 
the adaptive mesh. Thus, if this is the case, we add these neighbor cells to 
the adaptive mesh, nevertheless they are not evolved in time. Therefore, the 
memory ised is the sum of the cells of the adaptive mesh plus these neighbor 
cells. More details in 28, 27 . 

In order to control the L^-norm, Harten’s thresholding strategy is used, 
where 


= 0<i<L-l, 

|iZ| 


(7) 


and d = 2 is the space dimension and, in this two-dimensional case |ff| is the 
area of the domain. Therefore, in the Harten’s strategy, we use a smaller value 
of the parameter e in the coarser scales than in hues scales. For comparison, 
we shall also consider level independent threshold parameters: = e, for all 
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1. Herein, the multiresolution analysis corresponds to a prediction operator 


based on a third order polynomial interpolation on the cell-averages [27 


We recall that time integration is performed by a second order Runge-Kutta 
scheme. 


5. Numerical experiments 

We present here a 2D Riemann numerical experiment to illustrate the 
efficacy of our method compared to the traditional FV scheme. For the 2D 
Riemann initial condition we have used the values of the MHD variables 
presented in Table HJ The computational domain is [—1,1] x [—1,1] and 
Neumann boundary conditions have been applied. This example is proposed 
in j^, except for the boundary condition. 

We have also chosen 7 = 5/3, the final time of computations t = 0.1 and 
t = 0.25, the CFL parameter Ccfl = 0.3 and (^jch = 0.18. We have tested 
= e = 0.010,0.008,0.005 and Equation [7] with e° = 0.05,0.03,0.01. 

The reference GLM-MHD FV code used in this work has been developed 
in language, inspired by the Fortran code developed by j^, including 
an upgrade and new features for the implementation of the numerical flux 
HELD. The GLM-MHD MR code developed in [l^ is based on the hydrody¬ 
namics MR Garmen code developed in ^7, ^.The implementation has been 
optimized improving the momory allocation and unrolling the for-loops for 
the allocation of the variables. The GPU is improved about a factor 4for the 
test case studied here with L = 8 adaptive scales and e° = 0.01. 

For the numerical error analysis we have used a reference solution com¬ 
puted with a GLM-MHD FV scheme with L = 11 scales using the same 
numerical scheme in space, implemented in the AMROG code |8| which is 
parallelized. We computed the Li-error for the density solution [L\{p)). The 
GPU time for the MHD-FV reference is obtained with another code that is 
not parallel. 

The reference solution and numerical MR solutions for = 0.01 and 
L = 10 at f = 0.1 are presented in Figs. [2] and [3l respectively. For a later 
time t = 0.25, the numerical MR solution with L = 9 is presented in Fig. [71 In 
the solutions, we can observe that the structures are not always aligned, e.g., 
we can see a structure that appears in the density but not in the ^/-component 
of magnetic held in the right part of the domain. In this region, the latter 
variable is almost constant. This is expected because in plasma processes 
the discontinuities may not necessarily occur at the same position for all 
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Table 1: Initial condition of the 2D Riemann problem. The domain is [—1,1] x [—1,1] 


with Neumann boundary conditions and 7 


5 

3' 





X 

> 0 





y 

< 0 



y > 0 


p 

PUx 

pUy 

PUz 

P 

PUx 

pUy 

puz 

1.0304 

1.5774 

- 1.0455 

- 0.1016 

0.9308 

1.4557 

- 0.4633 

0.0575 

E 

Bx 

By 

B, 

E 

Bx 

By 

Bz 

5.7813 

0.3501 

0.5078 

0.1576 

5.0838 

0.3501 

0.9830 

0.3050 




X 

< 0 





?/ < 0 



y > 0 


P 

PUx 

pUy 

PUz 

P 

PUx 

pUy 

pUz 

1.0000 

1.7500 

- 1.0000 

0.0000 

1.8887 

0.2334 

- 1.7422 

0.0733 

E 

Bx 

By 

Bz 

E 

Bx 

By 

Bz 

6.0000 

0.5642 

0.5078 

0.2539 

12.999 

0.5642 

0.9830 

0.4915 


quantities. The component and p (not shown here) have a similar behavior 
as p, and the component has a similar behavior as B^. These observations 
are expected and they increase the number of cells in the adaptive mesh in 
the MHD case compared to hydrodynamic case. Fig. |4| presents an example 
of the adaptive mesh with e° = 0.01 for the initial, intermediate and hnal 
computational time. We can observe that the adaptive meshes represent all 
the structures present in the solutions. 

Using the GLM-MHD with the mixed correction, the divergence of the 
magnetic field is not necessarily zero. However, this correction improves the 
convergence of the numerical solution of the MHD system to the expected 
physical solution, as discussed in j^. Fig. [5] presents V-B for the FV reference 
for L = 11 and two MR solutions for L = 10 with e° = 0.01 at time t = 0.1 
and e° = 0.05 at time t = 0.25. We observe that the maximum values of 
divergence are in the front transition regions, near the central part of the 
domain. 

To check the time evolution of the divergence of the magnetic held, we 
consider the quantity 

Bdw{t) :=max{|V-B| : {x,y) G [-1,1]^}, 
where V ■ B is again evaluated using centered hnite differences. Fig. [HI shows 
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the time evolution of -Bdiv(^) up to t = 0.1 for the FV reference solution 
with L = 11 (d) and three series of MR computations with L = 8,9 ,10 (a, 
b, c) considering the following threshold, values e = 0, 0.010, 0.008, 0.005 
and e° = 0.050, 0.030, 0.010. For the reference solution we observe a rapid 
decay of the initial value, around 37, during the first iterations, followed by 
a relaxation towards the value 3 which is reached at about 0.04. Afterwards, 
this value remains almost constant. For the MR computations we find that 
not only the initial but also the relaxation values of Rdiv(^) depend on the 
hnest level L, and hence on the mesh size. For larger values of L the di¬ 
vergence becomes larger but in all cases we find that after a certain time 
-Bdiv(^) becomes constant or oscillates around a mean value. Using Harten’s 
strategy with e° these oscillations almost disappear. In Fig. |H] we consider 
the evolution of Bdiv{t) for longer times, up to f = 0.25, in MR cases with 
L = 9 for e = 0 and 0.005, and e° = 0.05. After t = 0.1 no oscillations can 
be observed for e = 0, while for both = 0.05 and e = 0.005 again some 
oscillations appear. 

One main conclusion in analyzing i?div(^) for the different cases is that 
no growth in time can be observed, thus the divergence error seems to be 
controlled by the divergence cleaning, as discussed in (2oj |. 

Considering the conservative quantities we compute the energy. 


S = 



(|up + |Bp) dxdy. 


and find the value 3.69 at the initial time. At time f = 0.1 we find for all 
FV solutions with L = 8,9 and 10 the value 3.48. For the MR computations 
we obtain 3.46, 3.47 and 3.48 for L = 8,9 and 10, respectively. These results 
are independent of the actual value of the threshold (ranging from 0.01 down 
to 0) and there is no signihcant influence if a fixed or level dependent value 
is used. This means that in all computations about 94% of the energy is 
conserved. At a later time, t = 0.25, we observe some decay, bnt still about 
86 % of the energy is conserved. 

The total magnetic helicity is also a conservative qnantity of the ideal 
MHD eqnations |3| and we consider its time rate of change, defined as. 


OH f f 

~ ® / / B ■ (u X B)dxdy. 


As shown in Fig. [HI right, the reference solution conserves perfectly the 
total magnetic helicity and dH dt yields values close to the machine precision. 


12 



For the three MR solutions there is an initial peak at about 4 ■ 10“^^ which 
immediately decays to near zero machine precision, and remains zero for 
e = 0. For the two others threshold values some intermittent spikes with 
amplitude below 2 • 10“^^ are observed. 



Figure 2: FV reference solution for the 2D Riemann problem using GLM-MHD with 
mixed correction. Shown are variables p, By, Uy and obtained at time t = 0.1 and 
L = 11. 

Table [2] presents a summary of the CPU time, memory compression, 
Dc and Tf(p) for all experiments at time t = 0.1. For = e = 0.005 
and e° = 0.05 the results are close, independent of the maximum level L. 
However, the case e = 0.005 has slightly better CPU time and memory 
compression with respect to L\{p). In these cases, for L = 10, the CPU time 
are 7 — 14% and the errors are approximately 10“^. As expected, the error 
increases for a scale-independent threshold = e with e being large, because 
it does not control well the error. However, as we decrease the value of e. 
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Figure 3: MR solutions with e° = 0.01 for the 2D Riemann problem using GLM-MHD 
with mixed correction. Shown are variables p, By, Uy and Uz obtained at time t = 0.1 and 
L = 10. 


the error becomes smaller. Thus, the choice of e is an important ingredient. 
We can observe that if we choose a sufficiently small e, both strategies will 
have similar behavior. However, we can optimize this process using Harten’s 
strategy, which corresponds to a level dependent e. 

In Table [3] we show the CPU time, memory compression, D^, and L\{p) 
for all experiments done at time t = 0.25. We present the simulations for 
= e = 0.005 and e° = 0.05. The results at t = 0.25 show that the 
MR approach does not introduce growing instabilities and it is possible to 
compute the solution for larger values of t. 
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t = 0 


t = 0.1 


t = 0.25 


1 

0.5 

0 

- 0.5 

-1 


-1 - 0.5 0 0.5 1 



-1 - 0.5 0 0.5 1 


Figure 4: Cell midpoints of the adaptive mesh L = 10 of the MR computation for the 2D 
Riemann problem using GLM~MHD with mixed correction at time t = 0 with 2.30% of 
the cells, at t = 0.1 with 26.65% and e° = 0.01; and at time t = 0.25 with 18.37% of cell 
and e° = 0.05. 



(a) (b) (c) 



Figure 5: Values of V • B for the 2D Riemann problem obtained with; (a) FV reference 
scheme using GLM-MHD with mixed correction for L = 11; and (b) MR scheme with 
e° = 0.01 using GLM-MHD with mixed correction for L = 10 at time t = 0.1; and (c) 
MR scheme with e° = 0.05 using GLM-MHD with mixed correction for L = 10 at time 
t = 0.25. Note that the values of this quantity are mesh-dependent. 
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(a) MR, L = 8 


(6) MR, L = 9 
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(c) MR, L = 10 


(d) FV, L = 11 




Time 


Time 


Figure 6: The quantity over time for the 2D Riemann problem, with; (a,b,c) 

GLM-MHD with mixed correction using the MR scheme with = e = 0.010, 0.008, 0.005 
and = 0, 0.05, 0.03, 0.01 for L = 8, 9, 10; (d) GLM-MHD with mixed correction using 
the FV scheme for the reference solution with L = 11. 
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Table 2: CPU time, memory, Dc, and density error L\{p) for the 2D Riemann prob¬ 
lem computed with MR scheme using GLM-MHD with mixed correction and either with 
constant or level dependent threshold for t = 0.1. 


L = 8 


= e 


MR 

€« 


FV 


0.01 

0.008 

0.005 

0.05 

0.03 

0.01 


CPU Time (%) 

22.74 

23.47 

24.55 

26.71 

27.80 

30.33 

100 

CPU Memory (%) 

44.18 

45.38 

47.70 

51.03 

53.12 

56.47 

100 

Dc (%) 

29.74 

30.67 

32.50 

34.94 

36.60 

39.28 

100 

Llip) -10-2 

3.680 

3.669 

3.657 

3.657 

3.652 

3.651 

3.640 

L = 9 


II 


MR 

eO 


FV 


0.01 

0.008 

0.005 

0.05 

0.03 

0.01 


CPU Time (%) 

13.63 

14.66 

15.91 

17.67 

19.00 

20.46 

100 

CPU Memory (%) 

27.03 

28.79 

31.24 

34.34 

36.01 

39.20 

100 

Dc (%) 

17.70 

18.97 

21.01 

23.51 

24.92 

27.42 

100 

Lfip) -10-2 

2086 

2.039 

1.981 

1.974 

1.958 

1.953 

1.9409 

o 

II 


= e 


MR 

€« 


FV 


0.01 

0.008 

0.005 

0.05 

0.03 

0.01 


CPU Time (%) 

7.73 

8.71 

9.85 

12.00 

13.03 

14.67 

100 

CPU Memory 

14.66 

16.02 

18.82 

22.40 

24.46 

27.48 

100 

Dc (%) 

9.25 

10.07 

12.01 

14.66 

1649 

19.25 

100 

Llip) -10-2 

1.090 

1.031 

0.932 

0.905 

0.895 

0.851 

0.841 


NOTE: The results are computed with second order Runge-Kutta for the MR scheme. 
The CPU time for the GLM-MHD FV method is 277 sec., 2326 sec. and 314 min., for 
L = 8, 9 and 10, at a Intel(R) Xeon(R) CPU E5620 2.40GHz, CPU 1596 MHz, cache size 
12288 KB and 4 cores. CPU time, memory and Dc performances are computed with the 
corresponding non-adaptive FV solution using L = 8, 9 and 10 scales on a uniform level. 
For the error, in all cases, we use a reference solution computed with a GLM-MHD FV 
scheme with L = 11 for the same numerical scheme, implemented in the AMROC code 

i- 
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Figure 7: MR solution for the 2D Riemann problem using GLM-MHD with mixed cor¬ 
rection for = 0.05. Shown are variables p, By, Uy and Uz obtained at time t = 0.25 and 


L = 9. 
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Figure 8: The quantities and time rate of change of magnetic helicity over time for 

the 2D Riemann problem, obtained with GLM-MHD with mixed correction MR scheme 
using = e = 0, 0.005 and e° = 0.05 for L = 9 and for reference solution. 


Table 3: CPU time, memory, Dc, and density error L\{p) for the 2D Riemann problem 
simulated with the MR scheme using GLM-MHD with mixed correction and with constant 
or level dependent threshold for t = 0.25 


L = 9 

MR 

= 0.005 eO = 0.05 

FV 

CPU Time (%) 

18.79 

22.61 

100 

Memory (%) 

38.12 

45.25 

100 

Dc (%) 

23.80 

29.03 

100 

Lfip) -10-2 

3.887 

3.826 

3.694 
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6. Conclusions and perspectives 

Starting from the ideal MHD equations completed with generalized La- 
grangian multipliers to control the incompressibility of the magnetic held, we 
have developed an adaptive multiresolution method in two space dimensions 
on a Cartesian mesh with local rehnement. The space discretization is based 
on hnite volumes with an HLLD numerical hux. For time integration an ex¬ 
plicit Runge-Kutta scheme has been applied. To introduce a locally rehned 
spatial mesh and also for local interpolation of the hux values Harten’s cell 
average multiresolution analysis has been used. 

To assess the efficiency and quality of this new adaptive scheme, we have 
considered a two-dimensional Riemann problem. We compared this numeri¬ 
cal solution with adaptive MR results for diherent threshold values and two 
strategies of varying resolution levels. The numerical results show that the 
divergence cleaning can indeed work successfully with adaptive space dis¬ 
cretizations. The MR method with constant thresholding exhibits better 
CPU time performance but worse precision when compared to the level de¬ 
pendent threshold. The only drawback with respect to the level dependent 
threshold computations is that the number of cells on the adaptive mesh is 
increased. We also observed that energy and time rate of change of mag¬ 
netic helicity, both conserved quantities in the ideal MHD equations, remain 
indeed approximately conserved in our adaptive MR computations. 

In future work we plan to complete the adaptive method with time adap¬ 
tivity using local and controlled time stepping and to perform thus fully 
adaptive simulations in three space dimensions. A second interesting direc¬ 
tion is to move to non-ideal MHD, taking into account resistive effects and 
hnite values of the huid viscosity to study the physics of reconnection of 
current sheets, especially in space physics applications. 
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Appendix HLLD Riemann Solver 

In the following solver, we consider the one-dimensional GLM-MHD equa¬ 
tions in their primitive form 


dE d 
dt dx 


^ dpu^ 
dt dx 




E -\- P ^ I Ux {^UxBx “t“ UyBy Bx 


d{pux 

dt 


A 

dx 

d{pUy 


pui+p + - Bl 


dt 

d{puz 


dt 


d 

(^pUxUy BxBy) 

d 

“f ^ {pUxUz BxBz^ 

dBx ^ 
dt dx 

3B d 

^ ~ ^xUy) 

dBz d , ^ ^ X 

^ ^ - BxUz) 

dip 2 dBx 
dt dx 


0 , (8a) 

0, (8b) 

0 , (8c) 

0 , (8d) 

0 , (8e) 

0 , (8f) 

0 , (8g) 

0 , (8h) 

( 8 i) 

cl 


Considering the MHD system described above, we can obtain the Jacobian 
matrix. From the structure of this matrix one can verify that the equations 
of Bx and ip can be decoupled from the remaining system and we can obtain 
the Jacobian matrix for the ID MHD system [^, p. 651-653]. The eigenvalues 
of this matrix are Ux, UxEcs, ± Ca and ± c/, where c^, c/ are the slow 
and fast magneto-acoustic waves and Ca is the Alfven wave. 

The Harten-Lax-van Leer-Discontinuities (HLLD) solver for MHD was 
hrstly developed by Miyoshi and Kusano and it can be considered as an 
extension of the Harten-Lax-van Leer (HLL) solver presented in [l^. The 
HLLD solver is based on four intermediary states and Q|j, 

divided by hve waves Sl, S^, Sm, S'^ and Sr, as illustrated in Fig. [9l These 
waves are related to the entropy, fast and Alfven waves. The HLLD nu¬ 
merical flux can resolve isolated discontinuities in the MHD system solution. 
This solver preserves positivity and it is more robust and efficient than the 
linearized solver, with an equally good resolution. 
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The states Q* and Q** for the GLM-MHD system are dehned as 


Q* = and Q^* = (p^*, P>r, B** 




with a denoting left (L) or right [K] states. In this approach, we compnte 
the nnmerical flnx for ip directly, then we consider p)* = pj** = ip here in the 
intermediary vector states, recalling that the HLLD is originally designed for 
MHD system, where the vector state Q has not the variable ip . 

The numerical flux function is given by 


/ 


Fhlld — ■* 


V 


Fl, if Sl > 0, 

FI, if < 0 < SI 
Ft, iiSl<0< Sm, 
Ft, if < 0 < St 

F|j, if < 0 < Sr, 
Fr, if Sr < 0. 


(9) 


The flux vectors Fz. = F(Qi), F^ = F(Q^) are exact, while F^, F|j 
are approximate fluxes at intermediary states Q^, Q^, and F^*, F^ are 
approximate fluxes at intermediary states Q^*, Q^. 

By the following process, we present the variables of the states Q* and 
Q**, allowing us to compute the HLLD flux in the intermediary states 

Fq = Fq, + Sa (Qq — Qa), /.|Q\ 

F- = F„ + Q- - - Sa) Q* - Qa, ^ ’ 

where a = R and L denote right and left, respectively. 

The following description of the HLLD flux is related to the x direction, 
considering B* = Bt = B^. In two-dimension, a similar expression can be 
obtained in the y direction, considering B* = B** = By. 

There are different possibilities to approximate the propagation speeds 
Sa', for instance, we use 


SL = mm{uL,UR)-max{cf^,Cft, Sr = max{uL,UR)+max{cf^,Cft, (11) 

where Ua are the plasma velocities, are the magnetic acoustic waves j^. 
The choice of 5^ is made to estimate the average normal velocity and it is 
given by 


Sm = 


{Sr — U^tlPR UxR — {Sr — U^tpL Uxl ~ PTr + PTl 
{Sr Uxt)PR {Sl UxiPPl 


( 12 ) 
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Figure 9: Schematic of the Riemann fan structure with four intermediate states used in 
the HLLD flux. Adapted from 221. 


The velocity is assumed to be constant over the Riemann fan, he., 




The total pressure pr = P + is kept constant, 

Pk = Pk = Pk = Ptu = Pk 


( 13 ) 


( 14 ) 


Under these conditions tangential and rotational discontinuities can be formed 
in the Riemann fan. 

From the choice oi Sm, the pressure can be written as 


Pt — 


(Sr - Ux^)pR PTl - [Sl - Uxl)pl PTr 
{Sr - Uxr)Pr - (Sr - Uxl)Pl 
PL Pr{Sr Rxr){Sl '^xl){^xr '^xl) 

{Sr - Uxr)pr - {Sl - Uxr)pl 


(15) 


Given Sm and p^, the states Q* = (p*,p*, 
are bordered by the states Qq, and they can be obtained from the jumps 
along Sa, where a = L or R represents the left or right state. Therefore, one 
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(16a) 


can derive the variables of the states Q* as 



6 

Pa 

— pa ^ 


~ '^Va 

< 

= 

Bl 

~ Bya 

BZ 

= B^^ 


a ^Xq 


Sa-Su' 


Sm ^Xo 


o ~ _ 

Pa{Sa-UxJ{Sa- Sm) - 

Pa{Sa - - Bl 


Pa{Sa - U^J{Sa “ Sm) “ 
Pa{Sa - - Bl 


‘ Pa{Sa - U^J{Sa “ Sm) “ B^' 
Consequently, we can compute 

{Sa - UxjEa - PTo,Ux^ + PtSm + Bx{\la ' B„ - U* ■ B*) 


Ei = 


Sa ~ Sm 


(16b) 

(16c) 

(16d) 

(16e) 

(17) 


During the computations some operations as 0/0 can appear when Sm = 
Ux^, Sa = Ux^ =t By^ = = 0 and B^ > 'ypa- In these cases, we have 

to replace u* = Uy , u* = u. ■, and B* = B* =0. 

^ Va yo! ’ 2 q: 5 ^q, Za 

Similarly, it is possible to obtain the equations related to the states 


Qa = iPa\Pa,<:,U;:,Ur 


B* 


TD-k-k 

' Vc. ’ 


BZ)- 


Due to the relation described by Eq. [131 starting with the jump condition 
of the continuity equation over an arbitrary value S', where Sl < S < Sm or 
Sm < S < Sr, we have 

pZ = pZ (18) 


The propagation velocities of the Alfven waves in the intermediary states are 
estimated by 


SI = Sm- 


\Bx 


V^’ 


S*r = Sm- 


\Bx 


\Zpr 


(19) 


Considering the jump conditions to the tangential components of the velocity 
and magnetic held over Sm, and if 7 ^ 0 , we can obtain the following 
relations 


u 

B 


•k-k 

Vl 

★★ 

VL 


** _ ** 

vr ~ y ’ 

D** — r?** 
Byu — By , 


ZD** ZD** _ ZD** 

B,^ = B^^ = B^ 
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( 20 a) 

( 20 b) 










If Bx = 0, it is impossible to calculate the remaining variables of the states 
Q**. Replacing Eqs. (TH] - [20] into the integral conservation laws over the 
Riemann fan, we can derive the variables 


u 


y 


u 


Z 


B 


y 


B 


2 


Kl + - Bl^)sigTl{B,) 

V^+ 

- ulJsignjB,) 

V/^ + 

+ \/pLpR{<a " ^z^)sign(^x) 


(21a) 

(21b) 

(21c) 

(21d) 


where sign(i? 3 ,) is 1 for B^ > 0, and —1 for B^ < 0. Consequently, the 
equation of the energy in Q** is given by 

K = -B; T VpI(U* ■ B; - u*‘ . B") sign(B,). (22) 

The same procedure is done for the y direction. 
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